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ABSTRACT 

We perform a series of high - resolution N-body simulations designed to 
examine the density profiles of dark matter halos. From 12 simulated halos 
ranging the mass of 2 x 10^^ ~ 5 x IO^^/i'^Mq (represented by ~ 1 million 
particles within the virial radius), we find a clear systematic correlation between 
the halo mass and the slope of the density profile at 1% of the virial radius, in 
addition to the variations of the slope among halos of the similar mass. More 
specifically, the slope is ~ —1.5, —1.3, and —1.1 for galaxy, group, and cluster 
mass halos, respectively. While we confirm the earlier simulation results that the 
inner slope is steeper than the universal profile originally proposed by Navarro, 
Frenk Sz White, this mass dependence is inconsistent with the several analytical 
arguments attempting to link the inner slope with the primordial index of the 
fluctuation spectrum. Thus we conclude that the dark matter density profiles, 
especially in the inner region, are not universal. 



Subject headings: galaxies: clusters: general - cosmology: miscellaneous 
methods: numerical 



1. Introduction 

Can one recover (some aspects of) the initial conditions of the universe from the 
distribution of galaxies at z ~ ? A conventional answer to this question is affirmative, 
provided that the effect of a spatial bias is well understood and/or if it does not significantly 
alter the interpretation of the observed distribution. This consensus underlies the 
tremendous effort in the past and at present to extract the cosmological implications from 
the existing and future galaxy redshift surveys. The two-point correlation function ^ (r) is a 
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good example supporting this idea; on large scales it is trivially related to the primordial 
spectrum of mass fluctuations, Pi{k). Furthermore the effective power-law index of the 
two-point correlation function on sufficiently small scales is related to the initial power-law 
index of Pi(fc) oc k""' as ^(r) oc r-3(3+"i)/(5+"i) (e.g. [Peebies(198(J]| ; [Suglnohara et al.(199l) 



Suto( 19931) . In other words, the initial conditions of the universe are imprinted in the 



behavior of galaxies on small scales (again apart from the effect of bias). This is why 
the phenomenological fitting formulae for the nonlinear power spectrum ( [Hamilton etj 
al.( 1991)1 , IJain, Mo, fc White (1995]1 , [Peacock fc Dodds(1996)| , [Ma (1998)|) turn out to be 
so successful. This fact, however, seems to be in conflict with the universal density profile 
proposed by Navarro, Frenk & White (1996,1997; hereafter NEW) for virialized dark matter 
halos. In their study, NFW selected halos which look to be virialized, and found that the 
density profiles universally obey the NFW form p(r) oc r~^(r -|- r^)"^. It is yet unclear to 
which degree their results are affected by their selection criterion which is not well-defined. 
In general, different halos should have experienced different merging histories depending 
on their environment and mass. Thus even if the halos do have a universal density profile 
statistically (i.e., after averaging over many realizations), it is also natural that individual 
halo profiles are intrinsically scattered around the universal profile ( |Jing(1999)| ). Definitely 
this is a matter of semantics to a certain extent; the most important finding of NFW is that 
such halo-to-halo variations are surprisingly small. 



A universal density profile was also reported by |Moore et al.(1999}| on the basis of 
high - resolution simulations of one cluster-mass halo and four galaxy-mass halos, and they 
claim that the density profile p(r) oc r^^'^ in the most inner region. In what follows, we 
will address the following quantitative and specific questions concerning the halo profile, 
especially its most inner region, using the high-resolution A^-body simulations; the inner 
slope of the halo profile is really described by p{r) oc or p(r) oc r~^'^ universally as NFW 
and [Moore et al.(1999]| claimed ? If not, does the slope vary among the different halos ? Is 
there any systematic correlation between the slope and the mass of halos ? 

In fact, some of the above questions have been partially addressed previously 
with different approaches and methodologies ([Fukushige fc Makino(1997)[ , [Evans[ 
fc Collet (199T][ , [Moore et al.(1998)[ , [Syer fc White(1998) 



Nusser fc 



Sheth(1999) , Jing(199?y , Avila-Reese et al.(1999)[) . In order to revisit those in a 
more systematic and unambiguous manner, we have developed a nested grid P^M N-body 
code designed to the current problem so as to ensure the required numerical resolution in 
the available computer resources. This enables us to simulate 12 realizations of halos in a 
low-density cold dark matter (LCDM) universe with (0.5 — 1) x 10^ particles in a range of 
mass lO^^^^^M©. 
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2. Simulation procedure 



As [t'ukushige fc Makino(1997]| and later |Moore et al.(199"8]| demonstrated, the inner 



profile of dark matter lialos is substantially affected by the mass resolution of simulations. 
To ensure the required resolution (at least comparable to theirs), we adopt the following 
two-step procedure. A detailed description of the implementation and resolution test will 
be presented elsewhere. 

First we select dark matter halos from our previous cosmological P^M N-body 
simulations with 256^ particles in a (lOO/i-^Mpc)^ cube ( |Jing fc Suto(1998)| , |Jing(1998)|) . 



To be specific, we use one simulation of the LCDM model of i^o = 0.3, Aq = 0.7, h = 0.7 and 
(Tg = 1.0 according to p:<-itayama fc Suto(1997)|. The mass of the individual particle in this 



simulation is 7 x 1O^M0. The candidate halo catalog is created using the friend-of-friend 
grouping algorithm with the bonding length of 0.2 times the mean particle separation. 

We choose twelve halos in total from the candidate catalog so that they have mass 
scales of clusters, groups, and galaxies (Table 0). Except for the mass range, the selection 
is random, but we had to exclude about 40% halos of galactic mass from the original 
candidates since they have a neighboring halo with a much larger mass. We use the 
multiple mass method to re-simulate them. To minimize the contamination of the coarse 
particles on the halo properties within the virial radius ai z = 0, rvir, we trace back the 
particles within 3rvir of each halo to their initial conditions at redshift z = 72. This is more 
conservative than that adopted in previous studies, and in fact turned out to be important 
for galactic mass halos. Note that we define r^ir such that the spherical overdensity inside 
is ~ IStt'^Qq^ ~ 110 times the critical density, Pcrit(-2 = 0). 

Then we regenerate the initial distribution in the cubic volume enclosing these halo 
particles with larger number of particles by adding shorter wavelength perturbation to the 
identical initial fluctuation of the cosmological simulation. Next we group fine particles into 
coarse particles (consisting of at most 8 fine particles) within the high-resolution region if 
they are not expected to enter the central halo region within 3rvir. As a result, there are 
typically 2.2 x 10^ simulation particles, ~ 1.5 x 10^ fine particles and ~ 0.7 x 10^ coarse 
particles for each halo. Finally about (0.5 — 1) x 10^ particles end up within r^i,. of each 
halo. Note that this number is significantly larger than those of NFW, and comparable 
to those of Ib'ukushige fc lVlakino( 1997)1 and |lVloore et al.(1998)| , |Moore et al.(1999)| . The 



contamination of the coarse particles, measured by the ratio of the mass of the coarse 
particles within the virial radius to the total virial mass, is small, about 10~^, 10"'^, and 
10"^ for cluster, group, and galactic halos respectively. 



We evolve the initial condition for the selected halo generated as above using a new 
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code developed specifically for the present purpose. The code implements the nested-grid 
refinement feature in the original P^M N-body code of [Jing fc Fang(199¥)] . Our code 
implements a constant gravitational softening length in comoving coordinates, and we 
change its value at z = 4, 3, 2, and 1 so that the proper softening length (about 3 times 
the Plummer softening length) becomes 0.004 r^r- Thus our simulations effectively employ 
the constant softening length in proper coordinates at z < 4. The first refinement is placed 
to include all fine particles, and the particle-particle (PP) short range force is added to 
compensate for the larger softening of the particle-mesh (PM) force. When the CPU time 
of the PP computation exceeds twice the PM calculation as the clustering develops, a 
second refinement is placed around the center of the halo with the physical size about 1/3 
of that of the first refinement. The mesh size is fixed to 360^ for the parent periodic mesh 
and for the two isolated refinements. The CPU time for each step is about 1.5 minutes at 
the beginning and increases to 5 minutes at the final epoch of the simulation on one vector 
processor of Fujitsu VPP300 (peak CPU speed of 1.6 GFLOPS). A typical run of 10^ time 
steps, which satisfies the stability criteria (Couchman et al. 1995), takes 700 CPU hours to 
complete. 



3. Results 

Figure |l| displays the snapshot of the twelve halos at 2; = 0. Clearly all the halos are far 
from spherically symmetric, and surrounded by many substructures and merging clumps. 
This is qualitatively similar to that found by [Moore et al.(199"8] , [Moore et al.( 199*9] for 



their high-resolution halos in the fig = 1 CDM model. The corresponding radial density 
profiles are plotted in Figure ^ The halo center is defined as the position of the particle 
which possesses the minimum potential among the particles within the sphere of radius 
around the center-of-mass of the fine particles. In spite of the existence of apparent 
sub-clumps (Fig. |I|), the spherically averaged profiles are quite smooth and similar to each 
other as first pointed out by NFW. The inner slope of the profiles, however, is generally 
steeper than the NFW value, —1, in agreement with the previous findings of [Fukushige fc| 
Makino (199711 and [Moore et al.(1998)| . We have fitted the profiles to p(r) oc r-^(r + r^)-3+/3 



with (3 = 1.5 (similar to that used by Moore et al.l999; the solid curves) and (3=1 (NFW 
form; the dotted curves) for 0.01r2oo r < r2oo, where r2oo is the radius within which the 
spherical overdensity is 200pcrit(-2 = 0). The resulting concentration parameter c, defined as 
the r2oo/rs, is plotted in the left panel of Figure |^. This is the most accurate determination 
of the concentration parameter for the LCDM model. There exists a significant scatter 
among c for similar mass ( [J ing(1999)[ ), and a clear systematic dependence on halo mass 
(NFW, [Moore et al.(TOTl ). 
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The most important result is that the density profiles of the 4 galactic halos are all 
well fitted hj (3 = 1.5, but those of the cluster halos are better fitted to the NFW form 
/3 = 1. This is in contrast with [Moore et al. (199^ who concluded that both galactic and 



cluster halos have the inner density profile p(r) cx r~^'^, despite that they considered one 
cluster-mass halo alone. In fact, our current samples can address this question in a more 
statistical manner. CLl has significant substructures, and the other three are nearly in 
equilibrium. Interestingly the density profiles of CL2 and CL3 are better fitted to the NFW 
form, and that of CL4 is in between the two forms. The density profiles of the group halos 
are in between the galactic and cluster halos, as expected. One is better fitted to the NFW 
form, whereas the other three follow the /? = 1.5 form. 

To examine this more quantitatively, we plot the inner slope fitted to a power-law for 
0.007 < r/r2oo < 0.02 as a function of the halo mass in the right panel of Figure |^. This 
figure indicates two important features; a significant scatter of the inner slope among the 
halos with similar masses and a clear systematic trend of the steeper profile for the smaller 
mass. For reference we plot the predictions for the slope, —3(3 + n)/{4 + n) by Poffman 



fc Shaham(1985) and —3(3 + n)/{5 + n) by Syer fc White(1998)| , using for n the effective 



power-law index n^g of the linear power spectrum at the corresponding mass scale (Table 
1). With a completely different methodology, [Nusser fc Sheth(1999)| argue that the slope 



of the density profile within ~ O.Olrvir is in between the above two values. On the basis 
of the slope - mass relation that we discovered, we disagree with their interpretation; for 
the galactic halos, the analytical predictions could be brought into agreement with our 
simulation only if the effective slope were —2, which is much larger than the actual value 
—2.5 on the scale. 

We would like to emphasize that our results are robust against the numerical resolution 
for the following reasons. Since we have used the same time steps and the same force 
softening length in terms of r2oo, the resolution effect, which is generally expected to make 
the inner slope of p(r) shallower, should influence the result of the galactic halos more 
than that of cluster halos. In fact this is opposite to what we found in the simulation. 
Furthermore our high-resolution simulation results agree very well with those of the 
lower-resolution cosmological simulations (open triangles) for the cluster halos on scales 
larger than their force softening length (short thin lines at the bottom of Figure §). We 
have also repeated the simulations of several halos employing 8 times less particles and 2 
times larger softening length, and made sure that the force softening length ~ 0.005r2oo 
(the vertical dashed lines of Figure 0) is a good indicator for the resolution limit. 
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4. Conclusion and Discussion 



In this Letter we have presented the results of the largest, systematic study on the dark 
matter density profiles. This is the first study which simulates a dozen of dark halos with 
about a million particles in a fiat low-density CDM universe. This enables us to address the 
profile of the halos with unprecedented accuracy and statistical reliability. While qualitative 
aspects of our results are not inconsistent with those reported by Moore et al. (1999), our 
larger sample of halos provides convincing evidence that the form of the density profiles is 
not universal; instead it depends on halo mass. Since mass and formation epoch are linked 
in hierarchical models, the mass dependence may refiect an underlying link to the age of 
the halo. Older galactic halos more closely follow the /3 = 1.5 form while younger cluster 
halos have shallower inner density profiles fitted better by the NFW form. Whether this 
difference represents secular evolution remains to be investigated in future experiments. 

Our results are not fully expected by the existing analytical work. Although the 
analytical work ( |Syer fc White( 1998)1 , [Nusser fc Sheth(1999)| , |Lokas(1999)| ) concluded 



that the inner profile should be steeper than —1, their interpretation and/or predicted 
mass-dependence are different from our numerical results. This implies that while their 
arguments may cover some parts of the physical effects, they do not fully account for the 
intrinsically complicated nonlinear dynamical evolution of non-spherical self-gravitating 
systems. 

We also note that the small-scale power which was missed in the original cosmological 
simulation has been added to the initial fluctuation of the halos. The fact that each halo 
has approximately the same number of particles means that more (smaller- scale) power has 
been added to the low mass halos than to the high mass ones. It is yet unclear how much 
effect this numerical systematics would have on the mass dependence of the inner slope 
found in this paper, and we will investigate this question in future work. 

In summary, the mass dependence of the inner proflle indicates the difficulty in 
understanding the halo density proflle from the cosmological initial conditions in a 
straightforward manner. Even if the density proflles of dark halos are not universal to the 
extent which NFW claimed, however, they deflnitely deserve further investigation from 
both numerical and analytical points of view. 
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Table 1. Properties of the simulated halos 



ID 






M 


a 




^vir 




GX 


1 


2, 


.30 


X 


10^2 


458,440 


0.269 


-2.47 


GX 


2 


5, 


.31 


X 


10^2 


840,244 


0.356 


-2.42 


GX 


3 


4. 


,21 


X 


10^2 


694,211 


0.330 


-2.44 


GX 


4 


5. 


,60 


X 


10^2 


1,029,895 


0.363 


-2.42 


GR 


1 


4, 


,66 


X 


10^3 


772,504 


0.735 


-2.31 


GR 


2 


4, 


,68 


X 


1013 


907,489 


0.736 


-2.31 


GR 


3 


5, 


,24 


X 


1013 


831,429 


0.764 


-2.30 


GR 


4 


5. 


,12 


X 


1013 


901,518 


0.758 


-2.30 


CL 


1 


4. 


,77 


X 


1014 


522,573 


1.59 


-2.12 


CL 


2 


3 


.36 


X 


1014 


839,901 


1.42 


-2.14 


CL 


3 


2, 


.89 


X 


1014 


664,240 


1.35 


-2.16 


CL 


4 


3 


.17 


X 


1014 


898,782 


1.39 


-2.15 



^Mass of the halo within its virial radius, in units 

of /i-IMq. 

'^Number of particles within its virial radius, 
'^the virial radius, in units of /i~iMpc. 

'^the effective slope of the linear power spectrum 
at the halo mass scale. 
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Fig. 1. — Snapshots of the simulated halos at ^ = 0. Left, middle and right panels display 
the halos of galaxy, group and cluster masses, respectively (see Table 1). The size of each 
panel corresponds to 2rvir of each halo. 
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Fig. 2. — Spherically-averaged radial density profiles of the simulated halos of galaxy (left), 
group {middle), and cluster (right) masses. The solid and dotted curves represent fits of 
/3 = 1.5 and (3=1 respectively (see text for the fit form). For reference, we also show 
p(r) oc r^^ and r~^'^ in dashed and solid lines. The vertical dashed lines indicate the 
force softening length which corresponds to our resolution limit. The open triangles in the 
right panel show the results for the corresponding halos in the lower-resolution cosmological 
simulation, and the long ticks at the bottom mark the force softening of the cosmological 
simulation. For the illustrative purpose, the values of the halo densities are multiplied by 1, 
10~^, 10~^, 10~^ from top to bottom in each panel. 
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Fig. 3. — Left panel: the concentration parameters for each halo for the Moore et al. 
(1999) form (filled circles) and for the NFW form (crosses). Numbers labeling each symbol 
correspond to the halo ID in Table 1. Right panel: Power-law index of the inner region 
(0.007 < r/r2oo < 0.02) as a function of the halo mass. The upper and lower dotted curves 
indicate the predictions of Hoffman & Shaham (1985) and Syer & White (1996), respectively. 



